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Abstract 

This paper discusses theoretical foundations oi 
quantitative image-based measurements tor extracting and 
reconstructing geometric, kinematic and dynamic 
properties of observed objects. New results are obtained 
by using a combination of methods in perspective 
geometry, differential geometry, radiometry, kinematics 
and dynamics. Specific topics include perspective 
projection transformation, perspective developable conical 
surface, perspective projection under surface constraint, 
perspective invariants, the point correspondence problem, 
motion fields of curves and surlaces, and motion equations 
of image intensity. The methods given in this paper are 
useful for determining morphology and motion fields of 
deformable bodies such as elastic bodies, viscoelastic 
mediums and fluids. 
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1. Introduction 

Image-based measurement techniques play an 
increasingly important role in virtually all natural sciences 
and engineering disciplines since they can provide 
tremendous information and knowledge about observed 
objects in a global, non-contact way with high temporal 
and spatial resolution. Specialists in photogrammetry. 
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computer vision, and other scientific and engineering 
disciplines have developed various methods that are best 
suitable to particular applications in their fields. In 
particular, both photogrammetrists and computer vision 
scientists have studied image-based techniques for many 
years to obtain metric and geometric information. The 
approaches developed by photogrammetrists are more 
mature and quantitative, which are recently extended to 
non-topographic applications [1]. By contrast, in order to 
deal with more complicated vision problems related to 
artificial intelligence, computer scientists tend to adopt 
more versatile mathematical approaches in perspective 
geometry, differential geometry and image algebra [2-5]. 
However, the approaches used by computer vision 
scientists are of qualitative nature in many cases and 
generally less accurate than those used in photogrammetry 
in metric measurements. Because the objectives of 
different disciplines are very different, there is a lack of 
sufficient interaction among specialists in various technical 
communities. Perhaps due to different notations, jargons 
and methodologies in these communities, it is difficult to 
transcend the different technical domains and see a unified 
scope of various image techniques. 

From a methodological standpoint, the approaches in 
photogrammetry and computer vision should be integrated 
into a universal theoretical framework. Furthermore, 
unlike computer vision scientists who mainly study rigid 
bodies, aerospace engineers and scientists often deal with 
complex morphology and motion fields of deformable 
bodies such as elastic bodies, viscoelastic mediums and 
fluids. It is highly desirable to formulate universal 
theoretical foundations for quantitative image-based 
measurements of morphology and motion fields of 
deformable bodies. In this paper, we will focus on the 
geometric, kinematic and radiometric aspects of image- 
based measurements. First, we will provide a unified 
treatment of the perspective projection transformation 
from the 3D object space to the 2D image plane and 
illustrate geometric connections among different 
formulations of the perspective projection transformation. 
Then, we will discusses some specific problems for 
recovering geometry and motion, such as projective 
developable conical surface, projection under surface 
constraint, reconstruction of motion field on a surlace and 
motion field of a 3D curve, the correspondence problem, 
and projective invariants. This is an area for combined 
application of approaches in perspective geometry, 
differential geometry, kinematics and dynamics. In the 
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radiometric aspect, we will discuss the fundamental 
relationship between the image intensity and radiance from 
an object. Based on this relation and imposed physical 
constraints, the motion equations of image intensity will be 
derived for typical physical processes such as moving 
Lamertian surface, emitting passive scalar transport, and 
transmitting passive scalar transport. These equations 
provide a rational way for reconstructing the geometric 
and kinematic properties of deformable bodies like fluids. 
In general, the geometric, kinematic and radiometric- 
approaches are closely coupled. 

2. Perspective Projection Transformation from 3D 
Space to 2D Image 

Image-based measurement techniques extract data 
from 2D images and then map them into the 3D object 
space. There is a perspective relationship between the 3D 
coordinates in the object space and the corresponding 2D 
coordinates in the image plane [ 1, 6-8]. Here, we discuss 
several formulations of the perspective projection 
transformation. Although these formulations are 
equivalent, one may be more convenient to use than others 
for a specific problem. The fundamental geometric 
problem in image-based measurements is to determine the 
object space coordinates X =( X \ X 2 , X ' ) T given the 
corresponding image (retinal) coordinates x =(x\x 2 ) 7 . 
Figure 1 illustrates the camera imaging process. The lens 
of the camera is modeled by a single point known as the 
perspective center (or the optical center), the location of 
which in the object space is X c =( X/, X ’ , X* ) r . 

Likewise, the orientation of the camera is characterized by 
three Euler orientation angles. The orientation angles and 
location of the perspective center are referred to as the 
exterior orientation parameters. The object space point, 
perspective center and image point lie along a straight line 
for a “perfect'' camera. This relationship is described by 
the collinearity equations, the fundamental equations of 
photogrammetry. On the other hand, the relationship 
between the perspective center and the image coordinate 
system is defined by the camera interior orientation 
parameters, namely, the camera principal distance c and 
the photogrammetric principal-point location 
x P = ( x' p . y * ) T . The principal distance e, which equals 

the camera focal length for a camera focused at infinity, is 
the perpendicular distance from the perspective center to 
the image plane, whereas the photogrammetric principal- 
point is where a perpendicular line from the perspective 
center intersects the image plane. Due to the lens 
distortion, however, perturbations to the imaging process 
lead to departures from collinearity that can be represented 
by the shifts Sx ! and Sx 2 of the image point from its 
“ideal" position on the image plane. The shifts Sx l and 


5x~ are modeled and characterized by a number of the 
lens distortion parameters. 

The image and object space coordinates of the points 
are related by the collinearity condition in which the image 
vector is aligned with the vector from the perspective 
center to the object space point 
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where M = / m y ] is the rotation matrix, X is a scaling 

factor. Algebraic manipulation of Eq. (2.1) yields the 
well-known collinearity equations (with the distortion 
terms Sx 1 and Sx 2 ) relating the point in the 3D object 


space to the corresponding point on the image plane, 

U_ C IL 

m s r (X-X e ) x J 

, o , m 2 r ( X - X ) ~X 2 
x -x,+ 8x‘ =-c— ^ — = -c— 

m 3 ( x - x ) j 


( 2 . 2 ) 


where the vectors m, =( ) T and 

tn 2 = ( m 2l ,m 22 ,m 2J ) T are the directional cosine vectors 
along the .v'-axis, .ir-axis in the image plane, respectively. 
The vector m 3 =( m 3f ,m 32t m 3J ) T is normal to the image 
plane, directing from the principal point to the optical 
center on the optical axis. As shown in Fig. I, the unit 
orthogonal vectors m,, m 2 , and m 3 constitute an object 
space coordinate frame at the optical center X c and 

X = ( X ,X ,X ) f are the projections of the object space 
position vector X - X c in this frame. The elements of the 
rotation matrix m tj ( i . j = j , 2, 3) are functions of the Euler 
orientation angles ( co t 0x ), 

m n = cos0 cos fC, rn /2 = sin co sin 0 cos k + cos CO sin K, 

m u = -cosco simp cosx + sinco sinx, m 2j = - cos0 sintc, 

tn 22 - - sin co sin 0 sin K + cos co cos x\ 

m 2 s - cos co sin 0 sin K + sincocostc, 

m 1t = sin0, = -sincocos0, m u =coscocos0. 


(2.3) 

The orientation angles ( co,0x ) are essentially the pitch, 

yaw, and roll angles of the camera in the established 
coordinate system in the object space. The rotational 
matrix M is an orthogonal matrix having the property of 

M 1 —M r or m? m j - S .. . The scaling factor 

X = -c/m/( X - X % ) is a ratio between the principal 

distance and the projected component of the object space 
position vector X - X c on the optical axis in -m 3 
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direction. When an object space point X is on the focal 


plane m/( X - X ) = 0 , the scaling factor becomes 
infinite, i.e., A = °o , which corresponds to the points at 
infinity on the image (retinal) plane. 

The terms Sx 1 and Sx 2 in Eq. (2.2) are the image 

coordinate shifts induced by the lens distortion. The lens 
distortion terms can be modeled by the sum of the radial 
distortion and decentering distortion [9-10] 

Sx l -8x[ + 8x ! (1 and Sx 2 =8x; + 8x] , (2.4) 

where, assuming that the optical axis of the lens is 
perpendicular to the image plane, we have 
Sx] = K f ( x 1 '-x' p )r + K 2 ( x 1 —x‘ j} )r 4 , 

Sx; = K,( x 2 ’-x; )r : + K : ( x 2 -x; }r 4 

Sx> = P, I r 2 + 2( x 1 -x' r ) 2 I + 2PJX 1 -x' if x 2 -x; ) , 

Sxj = P, [ r + 2( x 2 -x; ) 2 / + 2P l (x l -.v' )< x 2 -x 2 p ) 
r 2 =( x' +( x 2 -x; } 2 . (2.5) 

Here, K t and K 2 are the radial distortion parameters. P, and 
P 2 are the decentering distortion parameters, and .v' ’ and 
x 2 ' are the undistorted coordinates in image. When the 
lens distortion is small, the unknown undistorted 
coordinates can be approximated by the known distorted 
coordinates, i.e., x' —x 1 and x~ ~ x~ . For large lens 
distortion, an iterative procedure is employed to determine 
the appropriate undistorted coordinates to improve the 
accuracy of the estimate. The following iterative relations 
are used: ( x ' ’ /’ = .v ' and ( x' f = x~ . 

( X 1 ’ r ' = x 1 + Sx ' If x 1 ' ) k ,( X 2 ' / I and 

( x - ' ) M = x - + Sx 2 [( x 1 ’ / .< x 2 ’ ) k ] , where the 


superscripted iteration index k is k—0, /, 2 — . 

The col linearity equations Eq. (2.2) can be re-written 
in the homogenous coordinates in the image plane 
x„ =<x , .x 2 ,x 1 ) T =(x , .x 2 ,l) T 
Ax„ =AM<X-X C ) or x„=AP( X -X c ), (2.6) 


where P — [ p p ] — A 1 M and A — / a l} / is defined as 
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of Eq. (2.6) 


a xl = XrrijjfX’ -X‘ ), where the Einstein convention 

for summation is used. The matrix-form and tensor-form 
of the collinearity equations are sometimes convenient for 
mathematical manipulation. Another alternative form ol 
the collinearity equations in the homogenous coordinates is 
x h =AP h X h , (2.8) 


where X h = ( X 1 , X 2 , X ' . / / is the homogenous 
coordinates in the object space, and P h = A~‘M h and 
M h =(M - MX c ) are 3x4 matrices. Although Eqs. 
(2.6) and (2.8) are formally written as a linear relation 
between jt„ and X or i,, they are essentially non- 
linear because not only the lens distortion is a non-linear 
function of x . but also the scaling factor 

A = -dm , r ( X - X ) is not a constant in general. 

Nevertheless, because the lens distortion is usually small, 
its effect can be corrected by using an iterative scheme. 
Hence, Eqs. (2.6) and (2.8) can be treated as a quasi-linear 
system at each iteration. Without the lens distortion, the 
collinearity equations describe the ideal perspective 
projection. Eq. (2.8) is particularly suitable for utilizing 
useful results of classical perspective geometry to 
construct projective geometric invariants. 

Furthermore, Eq. (2.2) can be re-written as a form 
suitable to least-squares estimation for the object space 
coordinates X , 

W, T (X-X t ) = 0' (29) 

w/(X-X c i-o' 

where W, and W 2 are defined as 

(2 jQ) 

w 2 =(x 2 -x 2 p +Sx 2 )m 3 + cm 2 
As shown in Fig. 2. the vector W { is on the plane spanned 
by the orthogonal unit vectors m , and m ? , while W 2 is 
on a plane spanned by nt 2 and m 3 Geometrically 
speaking, W l r (X-X c ) = 0 and W 2 (X-X c ) = () 
describe two planes normal to W, and W 2 through the 

optical center. Thus, Eq. (2.9) defines an intersection of 
these two planes, which is a line through the optical center 
X c . For a given image point x = (x 1 , .v~ f , Eq. (2.9) is 
not sufficient to determine a point in the object space with 
the three unknown coordinates X =( X 1 , X ~ , X 4 ) T . 
Hence, extra equations associated with additional cameras 
and other geometrical constraints should be added for 
seeking a unique least squares solution of A . In contrast 
to Eq. (2.8), Eq. (2.9) does not include the scaling factor 
A . 

The collinearity equations Eq. (2.2) contain the 
camera parameters to be determined by geometric camera 
calibration. The parameter sets ( G),</>,!C,X J . X , X j , 
and (K r K 2 ,P l ,P 2 ) in Eq. (2.2) are the 

exterior orientation, interior orientation, and lens distortion 
parameters of a camera, respectively. Geometric camera 
calibration is a key problem in quantitative image-based 
measurements and a specific topic in both photogrammetry 
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and computer vision. Here we only briefly address this 
issue and readers can find the technical details of 
geometric camera calibration from references. In this 
paper, we generally assume that the camera is calibrated 
and a complete set of the orientation parameters and lens 
distortion parameters of the camera 
( **, X' , X ; . X;* . c. jr' ,x : ft , K , , K 2 ,P n P 2 ) is known. 

Analytical camera calibration techniques utilize the 
collinearity equations and distortion terms to determine 
these camera parameters [6-8]. Since Eq. (2.2) is non- 
linear, iterative methods of least squares estimation have 
been used as a standard technique for the solution of the 
collinearity equations in photogrammetry. However, 
direct recovery of the interior orientation parameters could 
be problematic and unstable since the normal-equation- 
matrix of the least squares problem is nearly singular. The 
singularity of the normal-equation-matrix mainly results 
from strong correlation between the exterior and interior 
orientation parameters. In order to reduce the correlation 
between these parameters and enhance the determinability 
of (c\x r ,y r ), Fraser |9, 11] suggested the use of multiple 

camera stations, varying image scales, different camera 
roll angles and a well-distributed target field in three 
dimensions. Nevertheless, the multiple-station, multiple- 
image method for camera calibration is not easy to use in 
many engineering and scientific applications like wind 
tunnel testing where optical access for cameras is limited 
and the positions of cameras are fixed. Abdel-Aziz and 
Karara [12] proposed a simple linear method for camera 
calibration. Direct Linear Transformation (DLT). 
Scientists in computer vision and robotics have developed 
various camera calibration schemes to achieve a fast 
calibration with an acceptable accuracy (a lower accuracy 
for a photogrammetric application). Tsai’s two-step 
method [13] is representative in computer vision, which 
uses a radial alignment constraint to obtain a linear least 
squares solution for a subset of the calibration parameters, 
whereas the rest of the parameters including the radial 
distortion parameter are estimated by an iterative scheme. 
By circumventing the singularity problem, Liu et al. [14] 
developed a robust optimization method for single-image, 
automatic camera calibration to determine the interior and 
exterior orientation parameters and lens distortion 
parameters plus the pixel spacing ratio. 

3. Projective Developable Conical Surface Containing 
3D Curve 

In this section, we introduce the concept of projective 
developable conical surface and show how r to reconstruct 
this surface containing a 3D curve from a single image. In 
principle, a 3D curve in the object space cannot be 
completely recovered from a single image since 
information in one dimension is lost in the imaging 
process. Nevertheless, using a calibrated camera, a 
projective conical developable surface on which a 3D 


curve lies can be reconstructed. When two calibrated 
cameras are used, the 3D curve can be uniquely 
determined as an intersection of two different projective 
conical developable surfaces. Furthermore, a 3D surface 
can be reconstructed as an envelope of a family of the 
projective developable conical surfaces obtained from 
images taken at different viewing angles. The motion Held 
of the 3D curve can be obtained from a time sequence of 
the curve. 

Generating Projective Developable Conical Surface 

Consider a 3D simple curve C in the object space, and 
its projection to the image plane and a plane P normal to 
the optical axis (parallel to the image plane), as shown Fig. 
3. The collinearity equations Eq. (2.6) are written as 

X-X C =A ‘Px h . (3.1) 

where P = P' = [ p.. J = M A = M T A . When the 

camera parameters and the scaling factor are constant and 
the lens distortion is fixed, differentiating Eq. (3. 1 ) yields 

dX = A 'P,: dx . (3.2) 

where dX = ( dX 1 , dX ’ . dX ' / . dx = (dx‘ ,dx 2 f . and 


(Pi, 



[ Pj, 


Pi! 


\ 


P22 

Pj2 


) 


A constraint imposed on Eq. (3.2) is m/dX-O, 

indicating that Eq. (3.2) actually describes the projection 
Cp of the 3D curve C on the plane P orthogonal to the 
optical axis direction or m } . This constraint is equivalent 
to the constancy condition of the scaling factor 
X = -c/m / ( X - X ^ ) since the differential 
dX - cm/ dX /[ m/( X - X c )] 2 shows 

dX —0<^>dX = 0. In fact, the constraint 

X~ - c / m/ ( X - X c ) = const . defines the plane 

orthogonal to the optical axis direction or m ? . As shown 

in Fig. 3, the projected curve C P on the plane P can be 
reconstructed from the image and then the developable 
conical surlacc D containing the 3D curve C can be 
generated. 

The arc length element of the projected curve C P on 
the plane P is 

dS CH = \dX I = X~ l \ Pv 1 1 ds , (3.3) 

where t — dx/ds and ds =1 dx I are the unit tangent vector 
and arc length element of the image of the 3D curve C in 
the image plane, respectively. Thus, the unit tangent 
vector of the projected curve C P on the plane P is 


dX 

dS r 

( p 


Pi2t 

i P <2 t I 


(3.4) 
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Note that the unit tangent vector T Cp is independent of the 

scaling factor A. The curvature vector of the projected 
curve Cp on the plane P can be obtained by differentiating 
Eq. (3.4) with respect to the arc length S Cp 


dT c 


K r = 
c ' - dS 


A 


( P<.k-T Cp 


d I Pj2 t I 
ds 


). (3.5) 


v, \P*tV 

where k -dt/ds = d 2 x/ds 2 is the curvature vector of the 


curve image in the image plane. The curvature vector k 
can be expressed as k =K c n , where /f and n ~k/ \ k \ 
are the curvature and the unit normal vector of the curve 
image in the image plane, respectively. Furthermore, we 


prove 

d I T a: 1 I _ (P.uk ) T ( Pnt ) 
ds \Pj2t\ 


(3.6) 


Hence, Eq. (3.5) becomes 
k A 




„ <Ps2nf(Ps2t), 

I Pa: n - T c = / • 


(3.7) 


I P i2 t\ 2 ' \Ps 2 t\ 

The curvature of the projected curve C P on the plane P is 
K Cp = K Cp • N Cp , where N Cp = K Cp / I K Cp I is the 

principal normal vector of the projected curve C/>. Thus, 
the ratio between the curvatures K Cp on the plane P and 


k on the image plane is 


x>„ 




I P a: t I 


1 /' in ^ I 1 


I P a: t I 


( 3 . 8 ) 


Clearly, Eq. (3.8) indicates that K Cp / K i is proportional to 
the scaling factor A . 

After the unit tangent vector T Cp is obtained from the 
image, the projected curve C P on the plane P is readily 
reconstructed by 

X c ,=X C f0+ )dS Cf . (3.9) 

J i) 


The initial position X Cp „ on the projected curve C P . in the 
object space is often chosen at the end point of the curve. 
Eq. (3.1) gives X Cr „-X t =r'Px he , where 
X M =(x' ir x-,l) T is the homogenous coordinates of the 
corresponding image point to X Crll . Substituting Eqs. 
(3.3) and (3.4) into Eq. (3.9) yields a ray vector directing 
from the optical center X c to a point X Cp on the 


projected curve C/> 

X Ct -X c =Ar'(Px h0 + f P.utds). (3.10) 
‘ J 0 

A family of the projective rays through the optica! center 
X c given by Eq. (3.10) generates a projective developable 
conical surface D that contains the 3D curve C. The 


tangent plane on the developable conical surlace D is 
given by 

(X-X e )*N D (s) = 0, (3.11) 

where N D ( s )- T Cp x( X Cp -X c )/\ T Cp x( X Cp - X c )\ 

is the unit normal vector to the tangent plane on the 
developable surface, which is independent of the scaling 
factor. Eq. (3.11) describes a single-parameter family of 
the tangent planes where the parameter is the arc length .v 
of the curve in the image plane. The projective conical 
developable surface, the envelope generated by the family 
of the tangent planes, is given by a system of Eq. (3.1 1) 
and Eq. (3.12) [15] 

(X-X c bdN D (s)/ds = Q. (3.12) 

Thus, the projective developable conical surface and 
associated geometric quantities such as the curvature, 
tangent vector and normal vector in the 3D object space 
can be obtained by using measured image quantities given 
the camera parameters. 

Reconstructing 3D curve and Surface 

From a single image, we are able to reconstruct the 
projective conical developable surface containing the 3D 
curve C rather than the 3D curve itself. Nevertheless, 
when two calibrated cameras are used, as shown in Fig. 4, 
the 3D curve C can be uniquely determined by intersecting 
the two projective developable conical surfaces associated 
with the different cameras. Interestingly, the developable 
conical surface intersection method for determining the 3D 
curve only requires knowing the correspondence of one 
distinguished point such as an end point of the curve. 

Furthermore, the developable conical surfaces can be 
used to reconstruct a 3D surface in the object space. As 
shown in Fig. 5, the developable conical surface 
containing the contour of the 3D surface can be 
constructed. Here the contour is a set of points on the 3D 
surface at which the surface normal is also the normal of 
the developable conical surface. When the camera is 
moved to a number of known positions through a 
rotational and translational transformation (rigid-body 
motion), a family of the developable conical surfaces can 
be obtained. The 3D surface is generated as an envelope 
of the family of the conical surfaces. Instead of moving 
the camera, the 3D surface can be rotated around a fixed 
axis such that a family of the conical surfaces can be 
obtained using a camera at a fixed position and viewing 
angle. From a computational viewpoint, this method may 
not be the most efficient since the intersection and 
envelope of the developable conical surfaces has to be 
determined. However, this method is to great extent 
immune from the ambiguous correspondence problem in 
stereovision. 

Recovering Motion Field of 3D Curve 

After two or more 3D curves in the object space at 
successive instants are reconstructed, we can estimate the 
motion field U( X ) of the 3D curve that is defined as 
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U(X ) = —. ( 3 . 13 ) 

dt 

The curve is given by X = X / S( t ),t ] , where / is time 
and S( t ) is the arc length of the curve in the object space. 
Measurements give the temporal and spatial difference 
between two curves at two successive instants t ( and t , 
(the time interval At = t 2 - t ! is small) 

A St X = XlS(t : )j 2 ] -X[ S( t f MJ. ( 3 . 14 ) 

Reconstruction of the motion field of the 3D curve from 
A st X is a non-trivial problem since the point 
correspondence between two sequential images is not 
known without using distinct targets on the curve 
especially for an elastic curve experiencing large and 
complicated deformation. 

The motion field of the curve is constrained by the 
underlying physical mechanisms behind the motion and 
deformation of the curve. In general, reconstructing the 
motion field is formulated as an optimization problem of 
the functional 

J[ U( X )] = 1 1 A s , X - U(X )At II -> min (3. 1 5) 
subject to relevant physical and geometric constraints 

G'[U(X )j=(), ( / = ) ( 3 . 16 ) 

and the suitable boundary conditions. Without the 
sufficient constraints, the solution to the optimization 
problem may not be unique. Also, the imposed physical 
constraints serve as a bridge connecting image-based 
measurements with the physical quantities in a specific 
problem being studied. 

In the simplest case in which the curve is rigid, the 
rigid-body motion field is expressed as 

U(X ) = U 0 +n„x(X-X 0 ) 9 ( 3 . 17 ) 

where U 0 and £2 (> are the constant translation velocity 
and angular velocity, respectively, and X 0 is the rotational 
center ol the curve. Because U 0 and £2 () together contain 
only six unknown constants, it is easier to solve the 
optimization problem. A slightly complicated case is that 
the curve is stretched in three fixed directions in addition 
to the constant translation and rotation. In this case, three 
stretching constants are added, and thus the total number 
of the unknowns in the optimization problem is nine. 
Next, we consider a highly deformable material line 
convected in an incompressible and irrotational flow. In 
this case, the physical constraints are the solenoidal and 
irrotational conditions [ 16] 

V*U( X ) = 0 and VxU( X ) = (). (3.18) 

A vortex-filament in an incompressible and irrotational 
How is an interesting example since the filament driven by 
not only mean flow, but also self-induction is no longer 
passive and the motion field is directly related to the 
geometric features of the filament. In this case, the 
induced motion velocity of the filament is proportional to 


the curvature K of the filament along the binomial 
direction vector B [ 17] 

U(X)ocfcB (3.19) 

Overall, the physical constraints for a specific application 
are necessary for recovering the correct motion field and 
associated physical properties of the 3D curve. 

4. Perspective Projection under Surface Constraint 

In general, mapping between a point in the 3D object 
space and the corresponding image point is not one-to-one. 
Nevertheless, as shown in Fig. 6, under a given surface 
constraint, a point on the surface has the one-to-one 
correspondence to the image point. In this section, we 
discuss the geometric relationship between the surface in 
the object space and the image plane. This topic is closely 
related to some applications in experimental fluid 
mechanics and aerodynamics such as reconstruction of 
complex flow topology from images of surface oil 
visualization and laser-sheet-induced fluorescence 
visualization. Consider a surface in the object space given 
by 

X< = F(X\X 2 j. (4.1) 

When Eq. (4.1) is imposed on Eq. (2.9) as a surface 
constraint, the perspective projection transformation Eq. 
(2.9) is reduced to 

( ^ / / ^ 2J )X “I" ( U’j VC, j W j j VI’,-! )X 

= w 2Jt W t T X t . W 2 t A 

w n X ’ + X - + IV, , F(X 1 . X 2 ) = W, T X . (4.2) 

where w t] ( i = 1,2 and j = 7,2,3 ) are the elements of the 
vectors W, =( w n , \v /2 , w IJt ) 7 and W 2 = ( w 2j , u\ : , \v, t f . 
For the given surface equation X J = F( X\ X 2 ) and the 
known camera parameters, the coordinates ( X 1 , X 2 ) T 
can be obtained from the image coordinates jc =(jc / , x 2 ) r 
by numerically solving Eq. (4.2). Thus, the coordinates 
X = ( X 1 , X - , X * ) T in the object space can be 
symbolically expressed as a function of the image 
coordinates x = (x* , x 2 f . that is, 

X =f s (x). (4.3) 

In fact, Eq. (4.3) is a parametric representation of the 
surface using the image coordinates x = (x 1 , .v* ) T as the 
parameters. Generally, the function f s ( x ) cannot be 
written as a closed-form solution except in some special 
cases such as a plane and a cylindrical surface. 

Differentiating Eq. (2.9), we have 

dW, T X +W, r dX =dW, T X c 
dW/x + W/dX = dW/X c . (4.4) 

When the lens distortion is fixed, dW, T =dx'm, ! and 
dW 2 T =dx 2 m/ hold. Then, substitution of 
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dX’=(dF/dX' )dX‘+(dF/dX 2 )dX 2 into Eq. (4.4) 
yields 


dX 


I \ 


dX 2 


= m/( X c -f s )Q 1 


dx' ' 
clx 


V 


where 


Q = 


w n + w , , dF / dX 1 w l2 + h / { dF / dX ~ 


(4.5) 


(4.6) 


vv, ; + u\, dF / dX 1 w\, + mv, dF / 3X ' 

Furthermore, the differential dX 1 can be expressed as a 
function of the image coordinates dx = (dx ! .dx 2 y 

dX ' = ( dF / dx 1 )d x‘ +(dF/dx 2 )dx 2 , (4.7) 

where 


dF 

dx" 


dF dX' dF dX 2 


ax' dx" ax- a.v 

Combining Eqs. (4.5) and (4.7). we have 
dX = m/( X c - f s )Qdx 


.(a = 1.2 


where 


Q 


Q 1 

dF dF j 

(-7T.-7T> /m 3 ( X c-Js , 

dx dx * ) 


(4.8) 


( 4 . 9 ) 


(4.10) 


Eq. (4.9) provides a fundamental relation between the 
differentials dX on the surface and dx on the image 
plane. The matrix Q is a function of the image 
coordinates, the camera parameters, and the geometric 
properties of the given surface. 

On the other hand, we notice 

dX =(dX /dx' )dx' +(dX / dx 2 )dx 2 . (4.1 1 ) 

From Eqs. (4.9) and (4.11), we obtain the following 
equality 

(dX/dx'.dX/dx 2 ) =m/(X ( -f s )Q. (4.12) 

The element dS of the arc length of a curve on the surface 
can be determined from Eqs. (4.11) and (4.12) from the 
image coordinates. We know 

dS : = I dX I = g a/) dx a dx t) , (4. 1 3) 

where 

ia 'p =L2) (4i4) 

is the so-called metric tensor in classical differential 
geometry [ 18). The summation convention is used in Eqs. 
(4.13) and (4.14). The quadratic differential form Eq. 
(4.13) is the first fundamental form of the surface in which 
the image coordinates are the parametric variables. In the 

case of the perspective projection transformation, X a/i 

may be properly named as the perspective metric tensor 
that is a function of the image coordinates, the camera 
parameters, and the properties of the given surface. 

The first fundamental form Eq. (4.13) allows us to 
measure the basic geometric quantities on the surface in 
the 3D object space from the image quantities. Consider a 


curve on the image plane given by a parametric form 
x(t) = (x'(t),x : (t)) T and the corresponding 3D curve on 
the surface X(t)= X( x(t))= f s ( x(l)) , where t is a 
parameter (e.g. time). The length of an arc bounded the 
points corresponding to the parametric values t = t„ and 


t =t, is 


S =\ U«p<dx“ /dt)(dx 


P /dt)J 1/2 dt . 


(4.15) 


The angle of two 3D curves at the intersecting point on the 
surface can be calculated based on the image quantities. 
Consider two image curves x(t) = (x l (t),x 2 (t)) 1 and 
x(t) = (x l+ (tl x 2+ (t)) T . The tangential vectors of the two 
3D curves on the surface are 
dX(x* ( t ), x 2 (t))/dt - dX/dx a dx a / dt and 

dX(x l+ (t), x 2 +(t))/dt = dX/dx a + dx a + / dt . Thus, the angle 
y of intersection is 

g aP (dx a /dt)(dx (i+ /dt) 

tOSY ~ Jy all <dx"/dl)<dx 11 /dt) j y aP (dx" + /dt)(dx /dr] 

(4.16) 

The area of a domain H on the surface can be expressed in 
the image coordinates 


A(H) 


if- " - 


(4.17) 


where U is the domain in the image (a^.av) plane 
corresponding to the domain H on the surface in the object 
space and g is the determinant g =1 g a/3 I . 

Example I: Plane 

The plane constraint is a simple, but very useful case 
in which the vector function f s (x), the matrices Q and 

Q can be explicitly expressed as a function of the known 
camera parameters and the measured image coordinates. 
Many aerodynamic flow structures are observed on a plane 
or a near-planar surface. Planar laser sheet How 
visualization is just a typical case of the plane constraint. 
In addition, a polyhedron consists of a number of the 
planar faces. Consider a plane in the object space 

X-* =a f X> +a 2 X : + a<. (4.18) 

This plane is defined by the vector a = (a r a 2 , a s f 
related to the normal vector of the plane. In this case, the 
matrix Q in Eq. (4.6) is 

+ vv /? a } w l2 + Wjj a 2 
^ w 2} + a f w 22 + w 2J a 2 , 

The function f s ( x) in Eq. (4.3) has a closed-form 
solution 


Q = 


( 4 . 19 ) 
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f six)-- 



where 


where 


/ = 


/ w l T x c 


a, 


(4.20) 


(4.21) 


(4.22) 


W 2 ‘X e - WlJ a s 
Now the matrix Q in Eq. (4. 10) is 

- ( Q~' 

Q = 

v ( a,, a, )Q , 

Example 2: Cylindrical Surface 

A cylindrical surface is another case where f s (x ) , Q 

and Q can be explicitly expressed. For the sake of 
convenience, a transformation from the Cartesian 
coordinate system to the cylindrical coordinate system is 
used, i.e., 

X =( X 1 , X \ X * ) T =( pcoscp , psitup , z) T , (4.23) 

where p is the radial coordinate, (p is the polar angle, and 
Z is the axial coordinate. The differential dX is 
f cosep - p sin (p 0 Y dp 3 
sin (p p cos (p 0 

\() 0 1 

For a cylindrical surface constraint p = p {) = const . , 
solving Eq. (2.9) for (p and z* we have f s (x) as a 
function of the image coordinates and camera parameters 


dX = 


A 


d(p 

dz 


(4.24) 


fs ( x > = ( Pa cos <p, p„ sill (p. Z f 


(4.25) 


where 


bjbi ± Jbjb; +b 2 -b;b 2 

cos (p — : : : 

bj+b; 

b 2 b , ± yjb;bj +/?/ — b 2 b] 


e= 


r w i:Pn cos (p - w,,p u sin (p Wj A 


(4.27) 


sm (p - 


W 2: p 0 cos (p - w 2 j p 0 sm (p w, ^ 

Another differential is dp = () . Note that the expressions 

of f s ( x j, Q and Q for a spherical surface can be also 
analytically derived, but they are so tedious that we do not 
present them here. 

5. Perspective Projection of Motion Field Constrained 
on Surface 

After discussing the geometric relationship between a 
surface in the object space and the image plane, we study 
kinematics under the surface constraint, that is, the 

perspective projection of a motion field on a surface. 

Consider a dynamical system 

^- = U(X), (5.1) 

dt 

where U( X ) = ( U r U 2 , U < ) T is a motion field in the 3D 

object space and t is time. A surface constraint imposed 
on the motion field Eq. (5. 1 ) is 

X* =F(X ! ,X 2 ). (5.2) 

Under this surface constraint, U( X j should be parallel to 
the surface, which obeys the orthogonality condition 

N x .U(X) = 0. (5.3) 

where N s ~ ( dF / dX 1 , dF / dX 2 , — / ) T is the normal 
vector of the surface. Under the surface constraint Eq. 
(5.2). Eq. (5.1 ) is effectively reduced to a 2D system 


(5.4) 


In fact, Eq. (5.4) describes an orthographic projection of 
the motion field Eq. (5. 1 ) onto the plane ( X 1 , X 2 ) . From 

Eq. (4.5), the dynamical system in the image plane, which 
is corresponding to Eq. (5.4), is 


d 

'X 1 ' 


'u J X' ,X 2 ,F(X’ ,X 2 )] ' 

dt 



U 2 l X‘.X : .F( X'.X : )] ; 


b] +/>; 


* 

d 

(x< > 

Q 

U,lfs<*)l ) 

p 0 sin (p — W t T X t . ) , 

u = 

dt 

A ) 

m/(X c -f s ) 

K U : ff s <x)l j 


h , =P„i ), 

b 2 = p„( u, : u\, - w„w tJ ) . 

b t = u\ ( W * X c - w lt W, T X c . 

There are two solutions for f s ( jr ), which are 
corresponding to two intersecting points between a 
perspective ray and the cylinder. For a non-transparent 
solid surface, a camera only sees one intersecting point at 
the surface facing the camera and hence / s ( x ) is one-to- 
one. The differentials in the cylindrical coordinate system 
are related to the image coordinate differentials by the 
following relation 

d( <p, Z) T - m/( X c - f s )Q-‘ dx , (4.26) 


(5.5) 


We call u = dx/dt = d/dt(x l y x 2 ) T the optic flow in the 
image plane. The optic How, a term first used in computer 
vision, is defined as the velocity field in the image plane 
that transforms one image into the next image in a 
sequence. If Eq. (4.2) gives a one-to-one topological 
mapping (homeomorphism): (x i , a ) h- > (X 1 , X 2 ) , the 

topological structure of the dynamical system Eq. (5.5) in 
the image plane is equivalent to that of Eq. (5.4) on the 
surface in the object space when Q has the full rank of 2 

and m/( X c - f s ) is not zero. Figure 6 illustrates this 
point. The problem is to recover two components of the 
motion field ( U r U 2 ) T using Eq. (5.5) from the measured 
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optic How u = dx/dt , while the third component Uj is 

readily obtained from the orthogonal condition Eq. (5.3). 

In the above analysis, we do not specify the motion 
field U( X ) , which could be a limiting viscous How field, 
an oil-film motion field driven by skin friction, or a 
particle motion field driven by a potential force (e.g. 
gravity and electromagnetic force). The physical 
constraints on U( X ), which are different in different 
cases, are necessary to reduce the number ol unknowns. 
For instance, an incompressible flow must obey the 
continuity equation 

V*V(X) = 0, (5.6) 


where V = (d /dX 1 t d/dX " ,d/dX ) T is the Laplace 

operator. Differentiating Eq. (5.3) with respect to X \ we 
have 


dU t 


dF dU, dF dU , 
- + ■ — " 


(5.7) 


ax-< dx 1 dx 3 ' dx 2 dx 3 * 

Substitution of Eq. (5.7) into Eq. (5.6) yields a constraint 
on (U ,, U , ) T for an incompressible flow field. 


3 dF a 
ax 7 + dx 1 ax' 




dF d 


dX 2 dX 2 dX 3 


- + 


U,=(). 


(5.8) 

In general, it is more difficult to directly obtain a 
global solution of Eq. (5.5) for the motion field 
(U U , f . Instead, we can seek a localized solution of 
Eq. (5.5) in a sufficiently small area. In a neighborhood of 
a point X (f , the motion field ( U ,, U 2 ) T can be expanded 
as a linear function of X 

U,( X ) = e m +( e „ , e i2 , e„ )( X -X„ ), (i = 1.2 ) (5.9) 

where e„, =U,( X 0 ) are the loeal velocity components 
and e :j = dU,( X 0 )/dX 1 (j = 1,2.3) are the local 

deformation components. Hence, the localized form of Eq. 
( 5 . 5 ) is written as 

Q 


Du=± 

dt 




m/(X c -f s ) 


( 


,+(*„. e l2 .e u )lf s (x)-f s (x„ )j 


(5.10) 


= 0 


\^ e 20 e 21. e 22 * e 2J )I f s( X ) f S ( X 0 H 

The unknowns e iH and e :j , can be determined by 
minimizing the norm II Du II , i.e., 

II Du II — » min . (5.11) 

At the final stage, the global motion field on the surface is 
reconstructed from the local motion fields. 

In an incompressible flow, the localized constraint Eq. 
(5.8) is 

dF dF 

e„ + e, , -^777 + e 22 + e 23 — - = 0 ■ 


dX 1 


dx- 


(5.12) 


Furthermore, for the irrotational motion field on a solid 
surface where the vorticity vanishes, i.e., 
a) = V x U( X ) = 0 , three constraints are 


U, 

U, 


d 2 F dF 

r + c,,— +(/. 


ax 'ax • 

a-F 


- + e, 


ax' 

dF 

ax' 


a-F a f 

t + e 22 — — - e 2 l = 0 . 


+ U. 


ax : ax 
a ; F 
ax -ax 


dX 
dF 

' +€:i ax- 


=o. 


ax 'ax' 

e 2l -e l2 =0. (5'3) 

Hence, for an incompressible, irrotational motion field, 
eight unknowns in Eq. (5.10) are reduced to four 
unknowns after these constraints are imposed. At the 
critical points, the velocity vanishes, i.e., 
e = UJ X„ ) = 0 . The local topological structures of the 
motion field at the critical points are determined by the 
deformation coefficients e n [19], 

The above method for calculating the local motion 
field is applicable to both discrete random particle patterns 
(e.g. particle image velocimetry (PIV) patterns) and 
continuous passive scalar patterns (e.g. laser-sheet-induced 
fluorescence patterns in fluids). When discrete particle 
patterns are so coarse that an individual particle can be 
tracked, the local optic flow u - dx/dt is the velocity ol 
the particle in the image plane [20-21 1. For dense discrete 
particle patterns, the local optic flow u = dx/dt can be 
obtained using PIV method to seek the maximum 
correlation between two particle patterns obtained at two 
consecutive instants. However, for continuous passive 
scalar patterns, recovering the local optic flow u = dx/dt 
w/h havp tn consider ihe oersnective 


nr\n_tri vifil 


c ini'i 1 


projection of the transport equations of passive scalar 
through a specific imaging process. Generally speaking, 
the perspective projection of physical processes will lead 
to motion equations of image intensity. The optic flow 
u - dx/dt is determined by solving the motion equation ol 
image intensity for a specific physical process given the 
suitable boundary conditions and constraints. Detailed 
discussion on motion equations of image intensity will be 
given in Section 12. 


6. The Correspondence Problem 

In Sections 4 and 5, three unknown coordinates in the 
object space are reduced to two when the surface 
constraint is imposed. Thus, the correspondence between 
the constrained surface and the image plane is one-to-one. 
In order to determine three unknown coordinates trom 
multiple views without any a priori constraint, however, 
we need to know the point correspondence between two or 
more images for the same physical point in the object 
space. This is the so-called point correspondence problem, 
one of the fundamental problems in 3D vision. Note that 
another correspondence problem is point correspondence 
in a time sequence ol images. Here we locus on the 
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stereoscopic correspondence of images rather than the 
temporal correspondence. 

Longuet-Higgins [22J gave a relation between the 
corresponding points in two images. Consider two 

cameras in which the unit vectors ( m lin) *m 2{n) , m 10i> ) 
constitute a local right-hand coordinate system whose 
origin is located at the perspective center X c(n) , where 

n- 1,2 is the index denoting the cameras 1 and 2. The 

three-dimensional coordinates X ini =( X{ n) , X ( 2 nr Xf n) ) T 

in the coordinate frames ( m llKr m 2w ,m Ma) ) are related by 

a tensor-form of the translation and rotation 

transformations 


*/;, = (6.1) 
where R = [ R afi j and T r = / Tf ] are the rotation matrix 

and translation vector, respectively. If the two cameras 
have the same principal distance and pixel spacing ratio, 
R and T r can be obtained by translating the origins X c(n) 

and rotating the vectors ( m Unr m 2(n) ,m Mn] ) (;; = /, 2) to 
match the two coordinates frames. Here R and T r are 

generally treated as the unknown matrix and vector. 

A new matrix Q is given by 

Q = RS or s at! = S /If 
where S is the skew-symmetric matrix 
'() Tf -T; N 
S=l-T/ 0 Tj 

t- - t; o 

Eq. (6.3) is written as a tensor notation 


( 6 . 2 ) 


(6.3) 


J 


S’ = f T° 


(6.4) 


where the permutation index e ^ = I, or - /, or 0 if 
( ) is an even, or odd permutation of ( J, 2.3 ) , or 
otherwise. From Eqs. (6. 1 M6.4), we know 


G*. *!n = KJ X? n -T r K )R„ U £ fl p„ r*7n 


)e M , a Trxf n =0 


(6.5) 


since R is orthogonal ( R up = S afi ) and is anti- 
symmetric in every pair of its subscripts. Note that 
X (n) =( X i , nr X ( 2 nr Xf ni ) r are the coordinates in the local 


frame ( m l(nr m 2(tl) whose origin is located at the 

perspective center. Thus, the collinearity equations Eq. 
(2.2) can be re-written as a simpler form. In the local 
coordinate frames ( ), without the lens 
distortion, the homogenous image coordinates 
I x h<nJ -( x ! nr x <n>'~ c ^ arc related to the object space 


coordinates X“ n by 

v « — — t ■ y / v 

X h<n> - C / A { „ 


(n = 1,2. a = 1.2,3) 


( 6 . 6 ) 


The image coordinates x“ nl are relative to the principal 
point in these local frames rather than the geometrical 
center of the image. Dividing Eq. (6.6) by Xf n Xf 2t /c 2 

yields the Longuet-Higgins equation for the image point 
correspondence 

xZ!,Q*4i, =0 - ( 67 ) 

Often, Q = lQ tifi ] is called the fundamental matrix that is 

related to the camera exterior orientation parameters. 
Given a number of the point correspondences between the 
two images (more than eight), the elements Q a(j can be 

determined by solving the following algebraic equations 
using a least-squares method 

( X hC) X h(l) h Qa0 = ^ ' ( / = L 2, ♦ • * ) (6.8) 

Longuet-Higgins’ original derivation of Eq. (6.7) is 
purely algebraic without giving a geometrical 
interpretation. In fact, the geometrical meaning of Eq. 
(6.7) is related to the epipolar lines in the images [2-3]. 
Given a point ( xj Ir x? n ) in the image 1, its epipolar line 

in the image 2 is a projection of the line connecting the 
object space point and the image point through the optical 
center in the camera 1 onto the image 2. The epipolar line 
in the image 2 is described by 

x t(2>Pam = 0 ' (6.9) 


where p ml) -Q af5 xf tlh are the coefficients of the epipolar 
line. Thus, the matrix Q maps the points in the image I to 
the epipolar lines in the image 2. in the same way, Eq. 
(6.7) also gives an epipolar line in the image 1 for a given 
point in the image 2. Hence, Eq. (6.7) serves as the 
epipolar constraint to reduce the number of unknowns in 
establishing the point correspondence. It is easily shown 
that when the lens distortion exists, the generalized 
epipolar constraint is 

( C, + : , )Qafl ( 4a, + 8xl n> = °- (6-10) 

The lens distortion terms are 

/ ^ x hno 1 = ( ^ x !nr^ x u,r^ ) 7 * Since the lens distortion 
terms in Eqs. (2.4) and (2.5) are non-linear, an epipolar 
line is a curve rather than a straight line. More point 
correspondences are required to solve Eq. (6.10) since 
there are additional unknowns associated with the lens 
distortion. 

The unknown fundamental matrix in the epipolar 
constraint is determined by using a number of point 
correspondences. Nevertheless, for two calibrated 
cameras, the image point correspondence can be directly 
established from the collinearity equations. The 
collinearity equations Eq. (2.9) for two cameras are written 
as 


W,J(X-X«, ) ) = 0 
W 2( /(X-X cl „,) = o' 


( 11 = 1.2 ) 


( 6 . 11 ) 
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Re-combination of Eq. (6.11) yields two sets of linear 
equations for X 


W lnm X=B lcom 

(6.12) 

W 2cum X = B lcom , 

(6.13) 


where the composite matrices and vectors are 



(w 7 1 

Vf 1(1) 


(w 7 ^ 

id) 

u = 

Icom 

w 7 

rv 2d) 

VV = 

' 2com 

W 7 

2d> 


w 7 

rr l(2\ 


w 7 

" 2(2) 


\ ) 


V I 





( \u 7 v ^ 

id) A cd) 

^ Icom 


' ^ 2com 

W T X 

2(1) A cd) 




W H2) X c , 2> ^ 


Eliminating X from Eqs. (6.12) and (6.13), we have a 
relation between the image coordinates (x l ur xf n ) in the 

image 1 and (xj 2r xf 2l ) in the image 2 

(*( Xj j ) i X~j ^ i ) ^ i com ^ 2 com ^ 2com ^ icom ^ 

(6.15) 

For a point (xl rr x; n ) in the image I, the corresponding 
epipolar line in the image 2 is given by 

\\G(x! !r xf n ;x! 2r x^)\\ = 0. (6.16) 

The Longuet-Higgins equation indicates that a point in 
the image 1 corresponds to the epipolar line on the image 2 
and vice versa. Therefore, the point correspondence is not 
uniquely established between a pair of images since given 
an image point there is only one equation for 

two unknowns (xj 2r x^)- In order to establish the point 

correspondence among images, we need at least four 
cameras (or lour images). For four cameras or images, the 
Longuet-Higgins equations are 

Qa,u H , 4» = 0 ,(i = 1, 2,3.4 ,j = l. 2,3.4 ) (6. 1 7) 

If the fundamental matrices Q„p u . jt are determined by 
calibration, for a given point (x[ n ,xf n ) in the image 1, we 
have a system of six algebraic equations for six unknowns 


( X (I) ' X(h ’ X(2i* 

*ai' x oo* 

[ U) ) 



X l,d> Qoc(U (-2 ) 

=°- 

X hd> Qap, 1-3 i X hU) 

= 0, 


X h{ 2) Qafit 2-3 > 

4s, =o. 

v« O X fi 

Xfild ( 1-4 f ' K h(4) 

= 0, 


X h(2) Qa/)( 2-4 ) 

5:^ 

■t. 

II 

X h(M Qapl 3-4 ) X h(4) 

= 0. 

(6.18) 


When the four cameras are suitably positioned. Eq. (6.18) 
is not singular and the solution of Eq. (6.18) for 
(x^.x^.x^.x^.x^.x^) can be obtained using an 

iterative method. In general, there are multiple solutions 
since three equations in Eq. (6.18) are quadratic. The 
correct solution has to be selected based on additional 


criteria. More than four cameras can be used to increase 
the redundancy for least square estimation. 


7. Composite Image Space and Object Space 

Eq. (6.12) gives a non-linear relation between the 
object space coordinates and X and the composite image 
coordinates x c „ m =(x' in ,x; l) .x l 0 _ t )' . As shown in Fig. 7. 
the local coordinate frame (m l<n ,m 2llr m 2fl) ) at the 
perspective center X c<l) on the image 1 can serve as a 
frame lor the composite image space in which 
x =(xh r xf lr xl 2f ) T are the coordinates along the unit 

vectors (m iw ,m 2ar m, (ll ). Note that the coordinate x,' 2l 
of the corresponding point in the image 2 is artificially 
assigned to the coordinate value in the axis rrt 3(l) in the 

composite image space. Mapping between the composite 
image space and the object space is one-to-one. 
Differentiating Eq. (6.12), we have 

W ICI , m dX+dW lcom X=dB lam . (7.1) 

Substitution of Eqs. (2.10), (6.12) and (6.14) into Eq. (7.1) 
yields a basic differential relation between the composite 
image space and object space (see Fig. 7) 
dX = HI x cnm )dx cnm or dX " = H ajS ( x am )dx ? om . (7.2) 


where 


H(x c , 


)=w,, 


'm M J(W lr : im B lc „ m -X r(n ) 0 O' 

0 m M J{W u ' nm B lam -X c(l) ) 0 

0 0 m l(2t T <W- c : m B lcom -X cl2) ) 

V ' 

(7.3) 


Consider a 3D curve in the object space. The arc 
length dS of the curve in the object space is expressed in 
the composite image coordinates, i.e. 

dS 2 = dX"dX a = J ap dx“ nm dxi ,„ , (7.4) 


where J aP = H pn H pp . Introducing the arc length 

ds =(dx“, m dx", m in the composite image space, we 


obtain a relation between dS and ds 

dS = Ux cnm )ds. (7.5) 

The length scale factor L( x cum ) is 

Ux a , m ) = (J aP C m t'i m ) ,n - , (7.6) 


where =dx"„„/ds is the unit tangent vector t com of 

the corresponding curve in the composite image space. 
Using Eq. (7.5), we are able to express the unit tangent 
vector T of the curve in the object space in the composite 
image space coordinates and tangent vector, i.e., 

T" = dX u /dS = L~‘ H a/) tH„ . (7.7) 

The principal normal vector K ot the curve in the 
object space is 
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K" = dT"/dS = L 


> H^L-' ) o 


a.vi 




+ L-'H n „ki 


all* ami 


(7.8) 

where A" m = dt“ im /ds = d 2 t" om / ds 2 is the principal 
normal vector of the corresponding curve in the composite 
image space. In the derivation of Eq. (7.8), the relation 
d / ds = t l * om d / dx ( * )t)) is used. The curvature of the curve 
in the object space is 

K <>bj =(K a K a ) ,/2 . (7.9) 

Eqs. (7.8) and (7.9) indicate that the curvature is not an 
invariant under the perspective projection transformation, 
which depends on not only k”,„. but also t“ om and the 

camera parameters. The unit principal normal vector jV is 
obtained by normalizing K 

N = k;1 K and N a = k;J„ K ° . (7.10) 

The unit hinormal vector B of the curve in the object space 
is 


B = T x N or B a = . (7.11) 

Thus, the torsion of the curve in the object space is 

r„, v =-N"dB" /dS 

(7 P) 

= -K-;J„K a L-‘e a/kl ( N»N a +T /I C m dN n /d.x?„ m ) 

In this stage, the geometric structures of the 3D curve such 
as the tangent, curvature and torsion are expressed as a 
function of the composite image space coordinates. In 
general, they are not differential invariants under the 
perspective projection transformation. In many 
applications, however, these geometric quantities are very 
useful since they are directly related to the physical 
properties associated with the curve. The useful physical 
properties can be extracted from them. For example, the 
motion of an isolated vortex filament (a good model for a 
tornado) is mainly determined by the curvature and torsion 
of the filament [17]. 

From Eq. (7.2), we can relate the motion field 
U (X ( X ) = dX a /dt in the object space with the motion 
field uj x com ) = dx* m /dt in the composite image space 

UJX) = H a/i (x com )u a . (7.13) 

The motion field U a ( X ) can be decomposed into two 
components 

UjX ) = dX a (S(t)j)/dt = dX a /dt + T a dS /dt . (7.14) 


The first term dX a /dt is the apparent velocity and the 
second is the deformation velocity along the curve. 
Similarly, uj x com ) has two components 

M ) = <KJs(t),t)/dt = dx^/dt + CJs/dt . (7.15) 
If the point correspondence of the curve at two successive 
instants is not known, Eq. (7. 13) cannot be directly utilized 
to calculate the motion field UJ X ) form image 

measurements. The deformation ds / dt in the composite 


image space cannot be determined from images without 
using any additional physical constraint. Thus, we have to 
look for a global method for recovering the motion Held 
that is briefly discussed in Section 3. 

8. Perspective Invariants of 3D Curve 

Construction of perspective algebraic and differential 
invariants for a 3D curve is difficult because the 
perspective projection transformation is non-linear. 
However, it is possible to construct semi-differential 
invariants in a special case of stereo image pair [231. The 
perspective invariants are useful since they can directly 
give certain geometric features of the curve from non- 
calibrated images. We use the perspective projection 
transformation for a pair of images 

X h(i) = \ i ) X h ,(/ = /, 2 ) ( 8 . 1 ) 

where x h(i) =(xj ir x f 2 ir I ) T is the homogenous image 
coordinates in a pair of images (/ = /,2), 
X h ~ ( X 1 . X “ , X \ 1 ) ! is the homogenous coordinates in 
the object space, and P h = (P^ m J ( n = I,2J< 
m = 1,2 J, 4 ) are a 3x4 matrix that only depends on the 
camera orientation parameters (see Section 2). In general, 
the scaling factors A (I , = -c U) / m 3r J ( X - X c(i) ) for the 

two images are not the same, which are related to the 
camera parameters and the position of a point in the object 
space. Here we consider a special but useful case in which 
the scaling factors in two images are equal, i.e., 

A<n = A,2 l =A- ( 8 . 2 ) 

The condition Eq. (8.2) implies 

c ,n=c, :i i m 3ll) = m JUj , m M J ( X c(2j - X cll) ) = 0 . (8.3) 
Eq. (8.3) indicates that the two images have a relative shift 
on the same plane normal to the vector m Ml) = m Mh . This 

means that two cameras are placed side by side and their 
optical axes are in parallel. This coplanar condition allows 
us to combine the collinearity equations Eq. (8.1) for the 
two images, which makes construction of perspective 
invariants possible. 

A relationship between the composite image space 
and the object space for a 3D curve is written as in the 


homogeneous coordinates 



* k c U JMS)) = MS)P hc , m 


(8.4) 

where X „,,,n \ r x! 2r l) T 

is 

the composite 

homogeneous coordinates in 

the 

image space. 

X h =( X 1 , X \ X \ / ) T is the homogenous coordinates in 


the object space, and P hcom is a composite matrix 
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hcom 


Phi nu 

Phi 1 H2 

Phi 1 )l< 

P \ 

r h{ / >14 

PfK 1 >21 

P/U i 122 

Phi l >23 

P 

r hi l >24 

Phi 2 Hi 

PfU 2 U2 

Phi 2 >13 

P 

r hi 2 >!4 

P/H 1 )M 

Phi 1 )32 

Phi 1 >33 

P 

‘hi t t34 j 


(8.5) 


The arc lengths s and S of the curves in Eq. (8.4) are used 
as a parameter of the curve in the composite image space 
(x‘ in , x; (r x{ 2i ) and the object space ( X \ X ‘ , X - ) , 

respectively. The function s = s( S ) is one-to-one. 

Brill et al. 1 23] has constructed projective invariants 
by differentiating Eq. (8.4) repeatedly with respect with S , 
arranging the results in matrix equations for several points 
on the curve, evaluating the determinants ol the matrix 
equations, and then eliminating all the factors related to the 
imaging parameters. At first, following the method 
developed by Brill et al. [23], we consider a number of the 
basic geometric structures. The curvatures in the 
composite image space and the object space are 1 1 8] 

K im = I x am \ = \x com x x am I and K„ bj = I X I = I X x X I , 

(8.6) 

where x mm = dx a)m / ds , x a>m = d • x nm / ds , 

X = dX / dS , and X =d : X/dS : are the derivatives 
with respect to the arc length. The torsions in the 
composite image space and the object space are, 
respectively. 


T im = I X rnm X rnm X com I / I X com 1“ 

^ hcom A hcom * hcom 


— 1 X X X h rom I / I X com I 


(8.7) 


and 

T obj = \X X *1/1* |- = -l *„ X„ 1/1*1-. ( 8 . 8 ) 

Eqs. (8.7) and (8.8) are expressed in the homogeneous 
coordinates x hcom = (*/,,, x‘, r x' (2) ,I f and 

X h -( X 1 ,X 2 ,X ! .1 ) T to facilitate the use of Eq. (8.4). 
In the object space, the unit tangent vector is T = * . the 
unit principal normal vector is N — X / , and the unit 

binormal vector is B - T x N = X x X / K ohj . We define 
the distance D :) from a point X, to the osculating plane to 
the curve at another point X j 
D tj = ( * j - X, ,). B j = I ( Xj - X, ) X j X, I / K ohj) 

= \X hJ X hj X hJ X hJ \ /K nh]j 

(8.9) 

where the subscripts ‘f and denote the quantities 
associated with the points *, and X j . The geometrical 

meaning of D :J is illustrated in Fig. 8. Similarly, in the 
composite image space, the distance d tJ irom a point 


x cami to the osculating plane to the curve at the point 

X 

** com. j 

d v = '(** 


-x 


hcom.i ) X hcom.j X hcom.j ^ im.j 


* hcom, j X hcom. j X hcom, j X 


I / K 
hcom, / im.j 


( 8 . 10 ) 


1 hcom.j 

= l 

In addition, we introduce the following geometric 
quantities 

i( 1 , 1 ,2, 3) I X X hcom, i X hcom, 2 X hcom, .f ^ ' 

i( 1 , 2, 2 , 3) = I X hcam J X hcom 2 X hcom.i X hcom, 3 ^ ’ 

m,r,2.3) = \X u X hJ *„, 2 **.,1, 

/(/, 2, 2\3) = \ * M X IU X„, 2 **..,!• (8-H) 

Differentiating Eq. (8.4) with respect to S , we obtain 


K hcom, i 


1 hcom. 


) 

ik \ 

( X>, A' A, > 

A i 


w ) 


*h,i 




X h,„m. I — ' PhamA X h.i ^ h.i ^ h.i ^ h.i ■ 


'L 

4„+L-2H2, 

+ A, + 2.VA, 

A,. 

( 8 . 12 ) 

where g,, = A, - A, s s~' • <?:, = -K ~ A, ■' s ' , and 

i = ds(S)/dS. From Eqs. (8.6)-(8.12), we have the 
following determinantal relations 


K~ T 

im.i on.i 


= A is* I P h 


K h j.i 1 

r 

' r*/>/ . i ’ 

(8.13) 

, 1 «■„/,. 

, D , ■ 

(8.14) 

^ Ph com 

1 /( IJ',2,3 ) , 

(8.15) 

1 Phcom 

1 1(1.2.2',3). 

(8.16) 


The subscripts T and denote the quantities associated 
with the points X i and X j in the object space and the 

corresponding points x comi and x com j in the composite 
image space. Re-arrangement of Eqs. (8. 13)-(8. 16) to 
eliminate , X 2 , s , and 1 P ) u . om I yields several semi- 
differential perspective invariants. 


(I) An invariant related to the torsions and the 
distances D n and d tJ is 

*im.l d 12 _ tpbjJ ^12 jyj 

T ,m.2 ^'21 T nhj.2 

For T nhj , = 0 , r ohi2 * 0 , D, 2 *0 f and D :I ±0. then 
r , ~0 . The zero-torsion point in the object space 

corresponds to the zero-torsion point in the composite 
image space. The condition D !2 * 0 , and l) 2} ^ 0 
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implies that the points X } and X 2 are not on the same 
osculating plane. 

(II) An invariant related to the curvatures, the 
distances D tJ and d t} . and the quantities 

HU\2J) and 1(L2 y 2\3) is 

* ia .2 d t2 r(U\2,3) r , D n I 2 (1J\2J) 

: ; = — — ^ . (8.18) 

K tm I d 2} r( 1,2,2' J) K oh!} D 2l 1-(L2,2\3) 

For K ou ,,2 = 0 . fC ()hj / * 0 , D r * 0 , and D 2} ± 0 , then 

x intm2 = 0 . This means that the zero-curvature point in the 

object space corresponds to the zero-curvature point in the 
composite image space. 

(III) An invariant related to the distances D tJ and d t; 

is 


d 2} d 4i _ P 2l l\ 
d 4i d 




(8.19) 


This result is analogous to the cross-ratio of the distances 
on a line, a classical perspective invariant in perspective 
geometry [2, 24]. 


9. Modeling of Imaging System 

Modeling of an imaging system is necessary for 
radiometric measurements. Figure 9 shows a radiation 
source at an infinitesimal area element dA y on the optical 

axis, having a distance R s from the optical center of the 

imaging system like a CCD camera 1 25 ) . The radiant 
energy (units: joule) from the area element integrated over 
a solid angle seen from dA s to a lens is 

dQ = dA s dt J U9,<j>)cos9dco , (9.1) 

where L(9 y 0) is the radiance (units: watt-m ^-sr' 1 ) of the 
radiative source at dA s , da)~ sin9d9d0 is the 
infinitesimal element of solid angle, 9 is the polar angle 
(measured from the surface normal), 0 is the azimuthal 
angle (measured between an arbitrary axis on the surface 
and the element of solid angle on the surface), and dt is a 
time interval. The number of photons collected by the lens 
is 

dn h „ iS =(hv)~ I T am dQ 
~(hvf I dA s dtT inm 

where hv is the energy of a single photon and T am is the 
transmittance of atmosphere air. Define 9 A as the angle 
between the optical axis and the line connecting dA s and 
the edge of the aperture. When 0 A is small, 
(sin9 A ) 2 ~ A„ / Rj is approximately the solid angle in 
which the radiative energy from dA s is collected by the 


j 


U 8,<p ) cos 6 dco 


(9.2) 


imaging system, where A„ is the imaging system aperture 
entrance area. Thus. Eq. (9.2) becomes 

= L p dA > d,T , < A, / Rj ) • (9.3) 

where L p is the average photon radiance over the 


collecting solid angle A n / R; ~(sind A f 

L p =(hv )~'(sin0 A r-J // 6,<p )cos0d(O . (9.4) 

Consequently, the number of photons reaching the image 
plane is 

= L p dA , dt ( A lt / Rj )T an „T ipl , (9.5) 

where T tpt is the transmittance of the optical system. The 

number of photons incident the detector element is simply 
proportional to a ratio between the detector element area 
A d and the image area dA } corresponding to dA s 

d ",,„ = d >'„„ A D /dA, . (9.6) 

Under the approximation of small angle 8 A « /, dA , is 
related to dA t by 


dA y / Rj = dA, / R \ , (9.7) 

where R : is the distance of the image plane from the 
optical center. Substituting Eqs. (9.5) and (9.7) to Eq. 
(9.6) and using the relations A„=nD : /4 and 
l/R, + l/R , — I / ft . we have 


dn 


(tt‘J 


n l -,. A„ dt T wll T >p , 

4 F 2 (l + M„ n ,f ' 


(9.8) 


where F = fl / D is the F-number and M = R, / R t is 

the optical magnification, D is the diameter of the 
aperture, and fl is the focal length. Thus, the total 

number of photons collected over an integration time t lST 
is 


= 


x L f , ^}) tfNT T am T opf 
4 F 2 U + M nnl ) 2 


(9.10) 


Since some of the variables in Eq. (9.10) depend on the 
frequency v of light, the number of photoelectrons 
generated in the solid-state detector over a frequency band 
[v r v 2 ] is 


V = f 

J V, 


, tt \ L (v )A t) t }NT T atm ( v )T ( v ) 

( V )— — ; — dV , 

4 F~( 1 + M Y 


(9.11) 

where R y (v) is the detector’s quantum efficiency (units: 


electrons/phonon). We separate the photon radiance L {> 
into the radiance magnitude L p independent of V and a 
shape function of the frequency spectrum f sp { v ) , i.e., 

L p =1 ~i\ r (v). (9.12) 

Thereiore, Eq. (9.1 1 ) becomes 
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v = /’.~v < 9I, > 

where p a „„ is a parameter describing the camera 
performance 

fJv)A„ ?/,vt T„J W )T„„,( v ) 




F 2 ( I + M r 


-dv . 


(9.14) 


After the camera is radiometrcally calibrated, the 


image intensity (gray level) is proportional to n f 
I(x) = c lm n lw . 


, i.e.. 


(9.15) 


The proportional constant c h „ is determined by 

calibration. The above analysis is made based on the 
assumption that the radiation source is on the optical line. 
In general, we have to take the off-axis effect into account 
[26-27]. Hence, a generalized form of Eq. (9. 1 5) is 

K x ) = c, p L cos 4 6 , (9.16) 

where 0 is the angle between the optical axis and light 
ray through the optical center. When the lens distortion is 
negligible, the angle 1 9 p can be expressed as a function of 

the image coordinates x , the principal point location x p 

and the principal distance c, i.e., 

6 v = arctcmi I x - x p \ / c). (9. 17) 

Grouping the terms in Eq. (9.16) that are only dependent 
of the image coordinates to the left-hand side, we get 

/( X ) OJ\ x-x p \) = c, m p ilim L p ( X ) , (9. 1 8) 

where the iunction describing the ofl-.ixis cited enn be 
approximated by 

OJ\x-x p \) = cos~‘0 p ~ l + 2\ x-x p l : / c- . (9.19) 

Assuming that the off-axis effect is corrected on the image 
plane, without loss of generality, we simply rewrite Eq. 
(9. 18) as __ 

l( x ) = c lm Pllim L p ( X ) . (9.20) 

In order to simplify the notations, we use replacements 
c„^c, mP , am and UX)^T p (X). Therefore, 

without loss of generality, Eq. (9.20) becomes 

l(x) = c m UX ), (9-21) 

where c m is a proportional constant related to the imaging 
system and U X ) should be understood as the spectrally 


averaged radiance. 


10. Typical Radiation Processes 

Surface Reflection 

Quantitative image-based measurements require the 
knowledge of the physical properties of radiation-matter 
interaction of the objects of interest. One of the important 
interactions is reflection on a surface. As shown in Fig. 
10, the incident radiance is generally a function of the 
incident direction ( ), i.e.. 


L t = L,( 6 r <p, ). (101) 

The reflection radiance L,(9 i .<p i ;B r .<p r ) is quantitatively 
characterized by the bidirectional reflectance distribution 
function (BRDF) [28] 

f r (0 i 4 l ;9,,<p r ) = dL r (9,,<t>,;9 r .<p, )/dE,(9,,tp i ). (10.2) 

where the infinitesimal incident irradiance dE t ( 9, ,0, ) 
over a solid angle element dco t is 

dE t f 9 t ,<p, ) = L i (d i ,0 i )cos9,dio : . (10.3) 

The BRDF has a unit of steradian The BRDF depends 
on the surface roughness distribution. Foe a perfectly 
diffuse surface or a Lambertian surface where the 
reflection radiance is isotropic, i.e.. L r = const . , the BRDF 
is /, = l/n . In this case, the reflection radiance is 

L r =( //7T )j L,(0, ,0, )cos9 l dco i . (10.4) 

( 0 , 

Furthermore, when the incident source of the irradiance 
E (i is collimated at a fixed incident direction ( 0 () ,<p o ) . the 
incident radiance is described by the Dirac-delta function 
,0, )=E„S(0 i -0 o )$(</>< -<p 0 )/ sin 0 n . (10.5) 

Thus, Eq. (10.4) becomes the Lambert's cosine law 

L r =( J / k )E 0 cos 0 () . (10.6) 

For a general surface, the BRDF can be derived based 
on either the wave equation for electromagnetic waves or 
geometrical optics. Using the method of Helmholtz- 
Kirchhoff integral, Beckmann and Spizzichino [291 have 
derived an expression for the mean power of 
electromagnetic wave scattered from a rough surface. 
Similar integral approaches were used by Icart & Arques 
[30) ar] d Wang [31]. Icart and Arques [30] derived an 
expression of the BRDF for multilayer materials, which 
was composed of specular, directional-diffuse (spread 
reflection), and uniform diffuse (Lambertian) components. 
From a viewpoint of geometrical optics, Torrance and 
Sparrow [32] gave a simpler expression for the BRDF. 
Beckmann-Spizzichino’s model and Torrance-Sparrow’s 
model were discussed by Nayar et al. 133 1 from a 
viewpoint of computer vision application. A 
bibliographical review on the BRDF was given by Asmail 
[34]. Scattering of electromagnetic waves from randomly 
rough surfaces is still an active research area covering a 
variety of theoretical and experimental studies [35], 

From a viewpoint of application, the empirical 
expressions for the scattered radiance from a rough surface 
are very useful due to their simplicity [36]. An empirical 
model for a single light source is 

L r ( X )= p a EJ X )+ p (i EJX )( N L s ) (10.7) 

+ Ph E,JX )p(R T V ) 

where the first, second and third terms are, respectively, 
the contributions from the ambient reflection, diffuse 
reflection, and specular reflection. In Eq. (10.7), p tl , p d , 
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and p s , are the empirical reflection coefficients for the 
ambient reflection, diffuse reflection, and specular 
reflection. As shown in Fig. 1 1, the vectors N , L s , R , 
and V are, respectively, the unit normal vector of a 
surface, the unit vector directing the light source from the 
surface, the unit main directional vector of the specular 
reflection, and the unit viewing vector. EJ X ) and 
E ( J X ) are the irradiances for the ambient environment 
and light sources, respectively. The function p( R T V ) is 
the directional distribution of the specular reflection, 
describing the spreading of scattered light. Phong [37] 
gave a power function p( R T V )-( R T V )\ In general, 
the main directional vector of the specular reflection, R , is 
a function of the incident direction of light - L s . Although 
there are theories for predicting R [32], it is not known for 
a general surface. The unknowns in Eq. (10.7), including 
R , the reflection coefficients and the parameters in 
p( R t V ), have to be determined by calibration. For 

multiple light sources, Eq. (10.7) includes superposition of 
the contributions from these light sources. 

Radiative Energy Transfer in Media 

When light travels in a medium, the radiance is 
affected by absorption, emission and scattering. The 
radiative energy obeys overall conservation of energy. 
The equation of radiative energy transfer can be derived 
based on a balance among absorption, emission and 
scattering, i.e., 

dL„ a C 

_ = S . VL, ; = S„ - /?„/.„ +-^J L n ( g, )0 n t S, .s )d(0, 

4n 

( 10 . 8 ) 

where s is the path vector, p n is the extinction 
coefficient, o srj is the scattering coefficient, 0 (s it s) is 
the scattering phase function, S 7 is a radiation source 
term, and the subscript rj denotes the frequency range. 
This transport equation has been used in radiative heat 
transfer [38] and radiative hydrodynamics [39]. Note that 
the terminology of the radiatve intensity (unit: 
watts/area/solid angle) used in literature of radiative heat 
transfer is just the radiance in radiometry. The solution 
techniques and the suitable boundary conditions have been 
discussed by Modest [38]. 

Luminescence 

Luminescence is an emission from molecules after 
they are excited by an excitation light with a suitable 
wavelength. Luminescent dyes, widely used as probe 
molecules in biological and medical applications [40], 
have been utilized for flow visualization and 
measurements. For example, based on oxygen quenching 
of luminescence, luminescent molecules immobilized in a 
polymer layer have been used for surface pressure and 
temperature measurements in aerodynamic testing. These 


new sensors are called as pressure- and temperature- 
sensitive paints (TSP and PSP). After luminescent 
molecules in PSP absorb the energy from the excitation 
light with a wavelength X h they emit luminescence with a 
longer wavelength X 2 due to the Stokes shift. Liu et al. 
[41 ] have analyzed luminescent radiation from a PSP layer 

and obtained the spectral luminescent radiance ( L ) 


L /2 = h 0( P , T )q 0 Es /2 U 2 )K 1 (p / /p)M(p ), (10.9) 

where 0( P,T ) is the luminescent quantum yield that 
depends on pressure (P) and temperature (7), Es / J/., ) is 

the luminescent emission spectrum, /? is the layer 
thickness, q () is the incident light flux, p = cos 0 is the 
cosine of the polar angle 0 , and the extinction coefficient 


P /{ = * s a Product of the molar absorptivity 8^ and 

luminescent molecule concentration c . The coefficient 
M represents the effects of reflection and scattering of the 
luminescent light at the wail. The term K l represents the 
combined effect of the optical filter, excitation light 
scattering, and direction of the incident excitation light. 

The luminescent irradiance E / n over a collecting solid 


angle Q is 


■I. 


L . cos 0 dQ 


= P, 1 h&( P. T)q 0 Es y Ja 2 ) K f < M >Q 


( 10 . 10 ) 


where < M > is the spectrally averaged quantity of M . 
Even though Liu’s analysis was focused on a thin PSP 
layer, calculation of luminescent radiance is generally 
valid for a luminescent volume where surface reflection is 
absent. The spectral luminescent radiance integrated over 
a volume V is expressed as 

k, = Esx A’ )K > V~'j 0 ( X )qJX ) X )dX . 


( 10 . 10 ) 

A similar analysis for the luminescent flux was given by 
Gaigalas et al. [42]. 


11. Reflection and Shape Recovery 

Reflection on a solid surface depends on the geometric 
properties of the surface. In principle, shape of the surface 
can be recovered from surface reflectance under certain 
conditions. Computer vision scientists have studied the so- 
called shape-from-shading problem for decades [43-44]. 
Here we give a general consideration that is particularly 
useful for more complex engineering structures. Figure I 1 
shows a surface element with the unit normal vector N . 
The incident polar angle 0. is the angle between the unit 
normal vector N and the unit vector L s directing the 
light source from the surface. The reflecting polar angle 
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6 r is the angle between the unit normal vector (V and the 
unit reflecting vector R . The azimuthal angle </> r is the 
angle between the projections of the vectors L s and R on 
the surface. Assuming that the reflecting vector R is on 
the plane spanned by L s and N , we have 

R = ct s ,N + a,L s . (HD 


The coefficients a N and a, are determined by solving the 
following equations 


cos 9 r — A r • R — a s + 

a L N • L s , 


R*R = a 2 N +2a N a { N 

• L s + a 1 ’ = I . 

(11.2) 

Eliminating a N from Eq. (11 

2) yields 


(7 + cos 2 9 i )a 2 L - 2cos 2 6 i a L 

+ cos 2 0 r -I=O. 

(11.3) 

There are two solutions for a 

i. 



cos 2 0, ± i[( I + COS 1 ' 9 t )(cos : 9 t --cos 9 r ) + I 
° l " I^cos r 6 { ~ 

(11.4) 

The reflecting polar angle 9 r is not necessarily equal to 
the incident angle 9 i especially at large incident angles 
due to the off-specular reflection phenomenon on a rough 
surface [32]. In general, 9 r >6 i insures that there is no 

imaginary solution for o l , which is also supported by 
experiment data. The condition 9 r >9 l indicates a L <0 , 
Thus, the appropriate solution for a N and a L are 

cos 2 9 i -^(I + cos 2 9, )(cos 2 6 , -cos 2 9 r ) + / 

Ql 1 + cos 2 9 { 

a N = cosO r -a L cos9 i . (4 1 -5) 

The reflecting polar angle 9 r can be expressed as a 
function of 9 t based on theories and experimental results. 
In a special, but very useful case 9 r = 9 i . Eq. (11.5) 
becomes 

cos 2 8 t — / _ ( (V • L s )~ - / 

1 I + cos 2 8j 1 + ( N • L s )~ 

2 cos 6 j _ 2 N • (116) 

" l + cos 2 9 l 1+(N*L , )-’■ 

Consider a surface X 1 = F(X'.X 2 ) illuminated by a 
single light source. The relation between the image 
intensity and reflection radiance from the surface is 
H x )=i\ n p a EJ X ) ( , 17) 

+ c E h ( X )l p d N . L s + p h p( R • V )J ' 

The relation between the image coordinates x and the 
object-space coordinates X is given by the collinearity 
equations Eq. (2.2). The unit normal vector N is 

N = (F x , . F x , - 1 f / + f 72 : + / . 0 1-8) 


where F x , =c)F/dX' and F x ,=dF/dX : . The unit 
vector L s directing the light source A",, from the surface 


L,=(X,-X)/ \X,-X\. (H.9) 

When the camera is sufficiently away from the object, the 
unit viewing vector V directing from surface to the 
camera is approximately 

V =-/«,, (1110) 

which is known for a photogrammetrically calibrated 
camera. The reflecting vector R is given by Eqs. (1 1.1), 
(1 1.5) and (1 1.6). Clearly, given an image intensity held 
l(x), Eq. (11.7) is a complicated non-linear first-order 
partial differential equation tor the surface 

X’ = F( X ' ,X 2 i . Thus, a numerical solution to Eq. 
(1 1.7) has to be sought with suitable boundary conditions 
and constraints. 

When the light source is away enough from the object 
relative to the size of the object, the incident irradiance 
EJ X ) and ambient irradiance EJ X ) can be 
considered to be homogenous on the surface ol the object, 
that is, EJ X ) = const, and EJ X ) = const. . In this 
case, the vector L s is also approximately homogenous and 
it becomes a constant vector. Thus, Eq. (11.7) is 


simplified to 

Ex) ( vvv P a E a ,| | ||^ 

+ c„ s EJ Pj N»L s +p , p(a N N »V + a L L s »V )] 

Eq. ( 1 1 . 1 1 ) is still complicated for analysis. Furthermore, 
at a Lambertian surface without the ambient illumination. 
Eq. (1 1.1 1) is simply 


/(Jt) = c E. p,, N • L t . 


( 11 . 12 ) 


In computer vision, a viewer-oriented coordinate system 
and orthographic projection are often used to further 
simplify the problem [45], The viewer-oriented 
coordinates ( X 1 ,X 2 ) in the object space are aligned with 


the image coordinates ( x 1 ,x~ ). The third viewer-oriented 

coordinate X s is in the direction of the viewing vector V . 
Eq. (11.12), known as the image irradiance equation in 
computer vision, has been extensively studied for shape- 
from-shading [43-44], For quantitative measurements, Eq. 
(11.12) can serve as the first-order approximation. 


12. Motion Equations of Image Intensity 

In this Section, we derive motion equations of image 
intensity from underlying physical principles. The motion 
equations of image intensity can be used tor recovering the 
optic flow and other physical properties from a time 
sequence of images ol continuous patterns. The temporal 
and spatial development of the image intensity depends on 
the radiation process that is characterized by the physical 
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parameters p = (/?,, p N ) T and the geometric 

parameters q =(q l ,q 2 r-- ,q M ) T , that is, 

Hx,t) = c^ U X, p,q,t ). (12.1) 

Differentiating Eq. (12.1) with time, we have the motion 
equation of image intensity 


a/ 

dt 


+ u • V . / = c 


dt * dt 


V ,L + ^-.V L 
r dt " 

( 12 . 2 ) 

where u = dx/dt is the optical How in the image plane, 
U = dX/dt is the motion field in the object space, and the 
gradient operators are defined as 
V x ^(d/dx\d/dx 2 
V x =(d/dX i t d/dX 2 ,d/dX J ) T , 
v /f =( a/d/v-.a/dp* /, 
v t/ =(a/a^,--*,a/a^, / . 

The first term in the right-hand side of Eq. (12.2) is the 
local temporal change of the radiance. The second term is 
the change induced by motion in a non-homogenous 
radiance field. The third and fourth terms are related to the 
changes of the physical and geometric parameters, 
respectively. Eq. (12.2) is a generic form of the motion 
equation of image intensity. However, the detailed 
structure of Eq. (12.2) depends on the specific physical 
process being studied. To determine the optical flow, 
Horn and Schunck [46] suggested the well-known 
brightness constraint equation dl /dr + u • V ^ / = 0 in 

computer vision. In fact, the brightness constraint 
equation is just an assumption that the image intensity 
remains invariant along a stream of images. Generally 
speaking, this assumption, which is not related to any 
physical process, does not hold exactly. In the following, 
we give the motion equations of image intensity for three 
typical cases. Similar results can be obtained for other 
physical processes. Determining the optic flow in the 
motion equation of image intensity is a constrained 
variational problem. 

Moving Lambertian Surface 

Consider a moving Lambertian surface illuminated by 
an incident irradiance field E,JX). Since the image 

intensity is l( x )- c sS!t E h p tl N • L s , the motion equation 

of image intensity for a Lambertian surface is 
dl 

a/ +M# ' 1 


= C,V Pd 


dN ^ 

<N.L s )(U.V x E h )+EJ . L s ) 

dt 


(12.3) 

The first term in the right-hand side of Eq. (12.3) is the 
change due to motion in the non-homogenous irradiance 
field. The term ( dN/dt)» L s represents the rate of change 


of the unit normal vector N of the surface projected in the 
illumination directional vector L s = ( L s , L s , , L s < ) T . We 

explore the connection of this term with the fundamental 
geometric quantities of the surface. The term (dN/dt)*L s 
is expanded as 


dN , dN r WI , 


(12.4) 


The surlace is described by a parametric equation 

* = (12.5) 

where and £ : are the parameters of the surface. The 
term L s *V v N can be expressed in £ 1 and 

L,.V x N= L sa - j^2 . . ( P = 1.2 , a = 1.2,3 ) (12.6) 

According to the formulae of Weingarten [ 1 8] 
dN _ dX 


df' 


= S aa b po 


d^ a ’ 


we obtain 


U.(L,.V X N ) = -l j} g aa b p 


dX 


(12.7) 


• U , (12.8) 


v* d g u 

where l p =L sa d^/d. X a , g™ are the contravariant 
metric tensor, and b fi(r are the coefficients of the second 

fundamental form of the surface. 

Emitting Passive Scalar Transport 

In a transport process of passive scalar such as 
fluorescent dye, scattering particles, and temperature in 
fluids, the radiance is assumed to be proportional to the 
density or concentration y/( X j ) of the scalar 

UX,t) = c ¥ y/(Xj), (12.9) 

where c v is a proportional constant. The density of the 
scalar y/( X,t ) obeys the transport equation 
dy/ _ dy/ 

~5T“~97 

where D y/ is the diffusion coefficient of the scalar. 

Differentiating Eq. (12.1) and using Eqs. (12.9) and 
( 1 2. 10), we have 

dl( x.t ) 


- + U .Vy/ = D w V 2 x y . 


( 12 . 10 ) 


dt 


-=c m c v D v V 2 x \f/ . 


(12.11) 


Furthermore, because of /( x,t )= c m c \/r( X ,t ), Eq. 
(12.11) becomes 

dl( x ,t ) 


dt 


~= D ^ : x I(X,t). 


(12.12) 


The Laplace operator V* can be expressed in the image 
coordinates x , i.e., 


jrj 2 d d 2 

X ~ hy dx r + hy ‘ dx a dx r 


(12.13) 
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where h x and h 
defined as 

K = 


\(X 


3 V 


(a = 1.2 ,y = 1.2 , P = 1,2,3 ) are 
dx r dx° 


and h = ■ 


(12.14) 


r dX^dX* “ >B dX * dX 11 ' 

For a photogrammetrically calibrated camera, h y and h n , 

are determined by the collinearity equations with the 
known camera parameters and the image coordinates when 
a surface constraint X * = F( X '.X 2 ) is imposed (see 
Section 4). Hence, the motion equation of image intensity 
for a passive scalar transport process is 

(12.15) 


a/ a/ _ n 

dt +Ua dx" * 


a/ . a 2 / 

r dx? + V dx a dx * 


The optical flow field u° = dx° /dt can be recovered 
from Eq. (12.15). In particular, using the orthographic 
projection x a = X ", we have 

/;„ = 0 and . (12.16) 


In this case. Eq. (12.15) is reduced to the standard 
diffusion equation [47] 


a/ a/ n ( 

— + «„-T — = D 
dt a dx " 


a 2 / 


(12.17) 


v dx a dx a 

Transmittant Passive Scalar Transport 

Here we derive the motion equation ot image intensity 
for transmittant passive scalar transport in a medium like 
fluids. When a light ray transmits through a bulk of 
passive scalar, the intensity of light is attenuated due to 
absorption and scattering, shown in Fig. 12. The radiance 
reaching a camera through the scalar medium is given by 

^L = S,WL = - PL. (12.18) 

ds 

where s is the path vector and (5 is the extinction 
coefficient. The solution of Eq. (12.18) gives the 
transmitted radiance 


L = L 0 exp( - I fids ) 




(12.19) 


Consider a bulk of the participating passive scalar confined 
by two virtual boundary surfaces T, and T, , as shown in 
Fig. 12. We assume that the camera is far enough away 
from the bulk of scalar such that the light path is almost 
parallel to the optical axis, i.e., s = -m , . In this case, it is 

convenient to use the object space coordinates X in the 
frame (m , ,m 2 ,tn , ), defined as 

x' = m, »(X -X c ) 

X 2 =m 2 .(X-X c ), (12.20) 

X * = m, .(X-X r ) 

where the unit vectors in , , m 2 , and in, are orthogonal, 
i.e., m pa m Ya = 8 yP . Under the above conditions, the 
transmitted radiance in Eq. (12.19) can be written as 


U X ,1 ) = L„ exp\ 


- f AX. 

Jr, 


t)dX 


( 12.21 


where the boundary surfaces are X = F ,( X ,X ,t ) and 
~x' = r,( H’dX'.t). The extinction coefficient is 
proportional to the concentration ipf X ,t ) of the scalar, 
i.e., 

P(X,t) = £ ¥ v(X,t), (12.22) 

where £ v is an absorption coefficient. The relationship 
between the image intensity and radiance is 

I( xj ) = c vvv L( X,t ), (12.23) 

where x=(x*.x 2 ) T is the image coordinates. 
Combination of Eqs. (12.21), (12.22) and (12.23) yields a 
basic relation between the image intensity and the 
concentration of the scalar 

• r, 

(12.24) 


/( x J ) — L n exp\ 


-e r f 


y/( X J )d X' 


Differentiating Eq. (12.24) with respect to lime, we 

have 
dl(xj) 


dt 


Hxd) 


(f 

, dr, 

A — =- - ^ 

. dr, " 

[Jr, dt 

u - dt 

' r > dt 


(12.25) 

Since y/( X ,t ) obeys the transport equation Eq. (12.10), 
the first term in the right-hand side is 

5‘V 




.dx‘ 

dX "3X " 


D. 


X -4 

Jr, dX° 


dX 


ax ax 

(a = 1,2,3) (12.26) 

The second equality in Eq. (12.26) can be easily proven. 
From Eq. (12.20), we know the differential relation 

d/dX" =m pa d/dX l> and then 

a 2 /ax"ax« = m /ia m rn d : /dxdx' (12 2?) 

-8 rP a 2 /d~x r d~x l> =a 2 /ax /, ax /i 

Integration by parts yields 

f a ~>„ d x* = — fVax' + B.T., 

Jr, dx ax ax ax Ji i 

(/} = 1,2. a = 1,2.3) (12.28) 

where the boundary terms B.T. are 
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B.T 




ox 


d£l 2 dv 


dX 


dx p 


or, 

ox p 


dw 

+ — — 

dr, dr, 

dy/ 

dr. or 


dx J 

r, dX P dX P 

ax' 

— ff — V 

OX OX 


4_ I//I 

vr, 

A _J 

2 r : . dr 

dy/ 

n ox p ox p 

n dx 

— fi h — j 

dx ax 

r . ax’ 


We consider that a bulk of the passive scalar is confined in 
a finite domain and the distribution of y/(X,t) rapidly 
decrease to zero outside the domain. This represents a 
typical case in many practical applications. Therefore, 
when the virtual boundary surfaces and /\ are large 
enough such that y/ and its derivatives at the surfaces 
approach to zero, i.e.. 


Hr, ~> 0 ' Hr. -» 0 • 


di// 


dx 


-»0, 


dy 


OX 


■0.(12.29) 


Since the boundary terms in Eqs. (12.25) and (12.28) 
vanish. Eq. (12.25) becomes 


dl( x,t ) 


= D v Ux.t) 


d 2 


ax^ax" 


s: 


if/ d X 


dt 

( 0 = 1.2) (12.30) 

Now we consider the transformation between the 
image coordinates x-(x‘,x 2 ) T and the object space 

coordinates X = ( X ' . X ' , X * ) J . The collinearity 
equations without the lens distortion are 


K —■< 


— xj! = —c X / X . (0 = 1,2) (12.31) 

Thus, from Eq. (12.30), the Laplace operator can be 
written as 

a- „ a-’ 


axVx" 


= A-' 


a x fi dx fi 


(0 = 1.2) (12.32) 


where A = — c/X is the scaling factor. Using Eqs. 
(12.24), (12.30) and (12.32), we obtain the motion 
equation of image intensity for transmittance images of 
passive scalar transport 


a/ 


a/ / 


d 2 I 


di a/ 


K dxW dx? dx p j 
(P = L2) ' (12.33) 

Note that a simple version of the motion equation of image 
intensity for transmittance flow images was given by 
Wilders et al. [48] based on the orthographic projection 
and other assumptions. 

13. Conclusions 

We study a number of theoretical problems in 
quantitative image-based measurements of geometric, 
kinematic and dynamic properties of observed objects 


(specifically deformable bodies). From a unified 
viewpoint, we discuss different formulations of the 
perspective projection transformation and their geometrical 
connection. These equivalent formulations of the 
perspective projection transformation are selectively used 
in this paper to study different geometric problems, 
depending on convenience of the formulation applied to a 
specific problem. The perspective developable conical 
surface containing a 3D curve is reconstructed from known 
image measurements of the curve. The developable 
conical surfaces can be used to reconstruct a 3D curve and 
a surface without solving the ambiguous correspondence 
problem in stereovision. Furthermore, the general 
methodology is proposed for reconstructing the motion 
field of a 3D curve from a time sequence of images. 

The perspective projection transformation under a 
surface constraint allows one-to-one mapping between the 
surface in the object space and the image plane. We 
explore the connection of the geometric structures and 
motion fields between the image plane and the surface in 
the object space. These issues are important in 
reconstructing the complex motion fields on a surface such 
as skin friction field on an aerodynamic body and passive 
scalar motion field illuminated by a laser sheet. Then, we 
consider the general point correspondence problem in 
multiple images. Longuet-Higgins relation for the point 
correspondence problem is generalized by taking the lens 
distortion effect into account. Generally, establishing the 
point correspondence requires at least four cameras or 
images. The concept of the composite image space is 
introduced. After the relationship between the composite 
image space and the object space is established under the 
coplanar condition, the perspective invariants of a 3D 
curve are constructed. These invariants allow us to 
directly know the geometric features of the curve such as 
torsion and curvature from images without calibrating the 
cameras. 

In the radiometric aspects, we discuss the relationship 
between the image intensity and the radiance received by a 
camera as well as typical radiation processes such as 
surface reflection, radiative energy transport through the 
participating mediums and luminescence. The motion 
equations of image intensity are derived for moving 
Lambertian surface, emitting passive scalar transport and 
transmittant passive scalar transport. These equations 
provide a rational foundation for recovering the optic 
flows and motion fields of deformable bodies (e.g. fluids) 
from a time sequence of images of continuous patterns. 
Future research will be focused on the development of the 
effective numerical techniques and algorithms and their 
implementation in various simulations and experiments. 
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Perspective center X { 











Figure 4. 3D curve as an intersection between two 
developable conical surfaces. 



Figure 5. Perspective developable conical surface 
containing a contour of a 3D surface. 


Figure 6. Mapping between 3D surface and image. 




Composite image space Object space 

Figure 7. Composite image space and object space. 



Figure 8. Geometric meaning of the distance D,,. 
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Figure 10. Incident light and reflection in a local 
coordinate system on a surface. 



Figure 1 1 . Vectors of incident, reflecting, and viewing 
directions. 






